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ABSTRACT 

O 

c$ ', A model that rapidly computes the secular evolution of a gravitating disk- 

planet system is developed. The disk is treated as a nested set of gravitating 
rings, with the rings'/planets' time-evolution being governed by the classical 
Laplace-Lagrange solution for secular evolution but modified to account for the 
disk's finite thickness h. The Lagrange planetary equations for this system yield 
a particular class of spiral wave solutions, usually denoted as apsidal density 
waves and nodal bending waves. There are two varieties of apsidal waves — long 
waves and short waves. Planets typically launch long density waves at the disk's 
nearer edge or else at a secular resonance in the disk, and these waves ultimately 
reflect downstream at a more distant disk edge or else at a Q-barrier in the disk, 
whereupon they return as short density waves. Planets also launch nodal bending 
waves, and these have the interesting property that they can stall in the disk, 
that is, their group velocity plummets to zero upon approaching a region in the 
disk that is too thick to support the continued propagation of bending waves. 

The rings model is used to compute the secular evolution of a Kuiper Belt 
having a variety of masses, and it is shown that the early massive Belt was 
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very susceptible to the propagation of low-amplitude apsidal and nodal waves 
launched by the giant planets. For instance, these waves typically excited orbits 
to e ~ sin 2 ~ 0.01 in a primordial Kuiper Belt of mass M^b ~ 30 Earth-masses. 
Although these orbital disturbances are quite small, the resulting fractional vari- 
ations in the disk's surface density due to the short density waves is usually 
large, typically of order unity. This epoch of apsidal and nodal wave propagation 
probably lasted throughout the Belt's first ~ 10 7 to ~ 5 x 10 s years, with the 
waves being shut off between the time when the large R > 100 km KBOs first 
formed and when the Belt was subsequently eroded and stirred up to its present 
configuration. 



1. Introduction 

The Kuiper Belt is a vast swarm of comets orbiting at the Solar System's outer edge. 
This Belt is comprised of debris that was left over from the epoch of planet formation, and 
this swarm's distribution of orbit elements preserves a record of events that had occurred 
when the Solar System was still quite young. A common goal of nearly all dynamical studies 
of the Kuiper Belt is to decipher this Kuiper Belt record. However this record is still open 
to some interpretation. 

The dots in Fig. 4 give the Kuiper Belt Objects (KBOs) eccentricities e and inclinations % 
versus their semimajor axes a. This Figure reveals the KBOs' three major dynamical classes: 
the Plutinos which inhabit Neptune's 3:2 resonance at a = 39.5 AU, the Main Belt KBOs 
which are the non-resonant KBOs orbiting between 40 < a < 48 AU, and the more distant 
Scattered KBOs that live in eccentric, nearly Neptune-crossing orbits. The Figure also shows 
that the Plutinos and the Scattered KBOs have inclinations that span < i < 30°, while the 
Main Belt KBOs appear to have a bimodal distribution of inclinations centered on i ~ 2° and 
i ~ 17° (Brown 2001). Note that accretion models show that these large ~ 100+ km KBOs 
must have first formed from much smaller planetesimal seeds that were initially in nearly 
circular and coplanar orbits having e and sini < 0.001 (Kenyon & Luu 1999). However 
gravitational self-stirring cannot account for the Kuiper Belt's current excited state, so one 
or more mechanisms must also have stirred up the Kuiper Belt since the time of formation. 

The orbits of the Scattered KBOs are perhaps the most easily understood. These objects 
likely have had one or more close encounters with Neptune which lofted these bodies into 
eccentric, inclined orbits (Duncan & Levison 1997). Repeated encounters with Neptune cause 
these objects' semimajor axes and eccentricities to evolve stochastically along the Neptune- 
crossing curve shown in Fig. 4, and most of these bodies are ultimately ejected or accreted by 
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the giant planets. However Neptune has numerous weak, high-order mean motion resonances 
that thread the Kuiper Belt, and these resonances permit some of these Scattered objects 
to diffuse to lower eccentricities. This allows a small percentage of the Scattered objects to 
persist over the age of the Solar System at eccentricities just below the Neptune-crossing 
curve seen Fig. 4 (Duncan & Levison 1997). This diffusion to lower eccentricities may also 
have been more pronounced had Neptune's orbit also migrated outwards. In particular, 
Gomes (2003) shows that during the epoch of planet migration, mean motion and secular 
resonances act as pathways that allow some scattered KBOs to descend irreversibly into 
lower-eccentricity orbits that are far from the Neptune-crossing curve and hence stable. 
Since Scattered KBOs have large inclinations of % > 10°, this process might also account for 
the Main Belt's bimodal inclination distribution, with the i ~ 2° component representing 
the Main Belt's native population of low-inclination KBOs, and the i ~ 17° being due to 
scattered invaders who were deposited in the Main Belt by Neptune. However this process 
is quite inefficient since only e ~ 0.1% of these Scattered KBOs manage to find orbits that 
are stable over a solar age (Gomes 2003). 

The possibility that Neptune's orbit had expanded outwards is also supported by the 
cluster of KBOs that inhabit Neptune's 3:2 resonance (see Figure 4). This may have oc- 
curred when Neptune had first formed and began to vigorously scatter the local planetesimal 
debris. This process can drive an exchange of angular momentum between the planets and 
the planetesimal disk (Fernandez & Ip 1984), so this episode of disk-clearing can result in a 
rearrangement of all of the giant planets' orbits over a timescale of ~ 10 7 years (Hahn & Mal- 
hotra 1999). Outward planet migration will also cause Neptune's mean-motion resonances 
to sweep out across the primordial Kuiper Belt, and these migrating resonances are quite 
effective at capturing KBOs and pumping up their eccentricities (Malhotra 1995). Models 
of this process show that if Neptune's orbit had in fact smoothly expanded some Aa ~ 7 
AU over a timescale longer than a few million years, then resonance capture would have 
deposited numerous KBOs in the 3:2 and the 2:1 resonances with eccentricities distributed 
over < e < 0.3 (Malhotra 1995; Chiang & Jordan 2002). Although numerous KBOs do 
indeed inhabit Neptune's 3:2 resonance, only a handful are known to live near the 2:1 reso- 
nance at a = 47.8 AU, and many of these bodies have eccentricities of e ~ 0.3 which puts 
them quite near the Neptune-crossing curve (see Fig. 4). Thus it is possible that some or 
all of the bodies orbiting near a = 47.8 AU might instead be members of the Scattered Belt. 
Thus if planet migration did indeed occur, then the apparent low abundance of KBOs at the 
2:1 resonance is a mystery that can only be partly due to the observational bias that selects 
against the discovery of lower eccentricity objects at the 2:1 (Jewitt, Luu, & Trujillo 1998). 

But of particular interest here is the Main Belt which preserves additional evidence for 
another mechanism having stirred up the Kuiper Belt. Again, if Neptune did indeed migrate 
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outwards Aa ~ 7 AU, then the entire Main Belt was swept by the advancing 2:1 resonance. 
Models of planet migration show that the efficiency of resonance capture is no more than 
~ 50% (Chiang & Jordan 2002), so the current members of the Main Belt evidently avoided 
permanent capture by slipping through the advancing 2:1 resonance. However this planet 
migration scenario is utterly unable to account for the high inclinations of < i < 30° 
observed in the Main Belt (see Fig. 4), since the N-body simulations show that the advancing 
2:1 resonance typically excites Main Belt inclinations by only a few degrees (Malhotra 1995; 
Chiang & Jordan 2002). Evidently an additional mechanism is also responsible for exciting 
the Main KBO Belt. It has been shown that sizable KBO excitation can occur if a recently- 
formed Neptune had been scattered outwards by Jupiter and/or Saturn into a Belt-crossing 
orbit (Thommes, Duncan, & Levison 2002). Another possible source of KBO excitation 
is the invasion of the Main Belt by the Scattered Belt (Gomes 2003). However this latter 
process is a very inefficient mechanism having e ~ 0.001; an alternate mechanism that is 
possibly more efficient will be explored below. 

It has also been suggested that secular resonance sweeping may have been responsible 
for exciting the Kuiper Belt (Nagasawa & Ida 2000). Note that the locations of the secular 
resonances are very sensitive to the Solar System's mass distribution. Consequently, the 
depletion of the solar nebula (which includes perhaps ~ 99% of the Solar System's initial 
mass content) could have driven these secular resonances across vast tracts of the Solar 
System. Indeed, the model by Nagasawa & Ida (2000) suggests that if the nebula was 
depleted on a timescale of t ~ 10 7 years, Main Belt inclinations of i ~ 20° can be excited 
due to the passage of the secular resonance that sweeps outwards to infinity as the 
nebula is depleted 1 . However this model is rather idealized in that it treats the Kuiper Belt 
as massless, which is a concern since a primordial Kuiper Belt having some mass is also 
susceptible to the propagation of very long-wavelength spiral waves that can be launched 
at a secular resonance in the disk (Ward & Hahn 1998, 2003). This issue is worth further 
examination since wave-action can alter the magnitude of resonant excitation considerably 
(Hahn, Ward, & Rettig 1995; Ward & Hahn 1998). But even if there is no resonance in 



1 It should be noted that the resulting inclination excitation is also sensitive to the tilt between the 
nebula midplane and the invariable plane. For instance Nagasawa & Ida (2000) place the nebula midplanc 
in the ecliptic, and this results in substantial excitation. But if the nebula midplane is instead placed in 
the invariable plane, which is tilted 1.6° from the ecliptic, then almost no excess excitation results (Hahn & 
Ward 2002). We also note that the nebula models of Nagasawa & Ida (2000) as well as Hahn & Ward (2002) 
both treat the gas disk as a rigid slab of gas. However a more realistic treatment would allow the nebula 
disk to flex and warp in response to the planets' secular perturbations. It is suspected that this additional 
degree of freedom will substantially alter the secular resonance sweeping; indeed, it can be argued that the 
^15 never did sweep across the Kuiper Belt on account of this flexure (E. Chiang and W. Ward, 2002, private 
communication), so perhaps secular resonance sweeping of the Kuiper Belt is actually a moot issue. 
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the disk, planets orbiting interior to a particle disk can still launch these spiral waves at 
the disk's inner edge (Ward & Hahn 1998). Preliminary results from a model of secular 
resonance sweeping also reveals that a Kuiper Belt having only a modest amount of mass 
is utterly awash in these waves once the nebula is depleted (Hahn & Ward 2002). However 
the purpose of the present study is to first characterize the properties of these waves in the 
simpler, post-nebula environment, and to explore their cosmogonic implications. Thus the 
following will consider a suite of models of the secular evolution of the outer Solar System 
for primordial Kuiper Belts having a variety of masses. 

Accretion models tell us that the primordial Kuiper Belt must have had a mass of 
Mkb ~ 30 M® in the 30 < a < 50 AU interval in order for Pluto and its cohort of KBOs to 
have formed and survive over the age of the Solar System (Kenyon & Luu 1999). A similar 
Kuiper Belt mass is also needed to drive Neptune's orbital migration of Aa ~ 7 AU (Hahn & 
Malhotra 1999). However the current mass is M KB ~ 0.2 M e (Jewitt, Luu, & Trujillo 1998), 
so the Kuiper Belt appears to have been eroded by a factor of ~ 150. This may be due to 
a dynamical erosion of the Belt by Neptune or possibly by other perturbers that may once 
have been roaming about the outer Solar System, as well as due to the collisional erosion 
that has since ground all of the smaller KBOs down to dust grains that are then removed 
from the Solar System by radiation forces (Kenyon & Luu 1999; Kenyon & Bromley 2001). 
The following will give results obtained from models of the secular evolution of the outer 
Solar System for Kuiper Belts having masses in the interval < M K b < 30 M®. 

Section 2 derives the so-called rings model that will be used to study the secular evo- 
lution of disk-planet systems; the reader uninterested in these details might skip ahead to 
Section 3 or 4. Since spiral density and bending waves appear prominently in the model 
results, their properties are examined in Section 3. Section 4 will describe the model's ap- 
plication to the primordial Kuiper Belt, and a summary of results is then given in Section 
5. 



2. The Secular Evolution of Disk-Planet Systems 

The following shall treat the disk as a collection of nested gravitating rings in orbit 
about the Sun. Their mutual perturbations will cause these rings' to slowly flex and tilt over 
time, and this evolution is governed by the Lagrange planetary equations. 
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2.1. the rings model 

Begin with the gravitational potential that a single perturbing ring of mass m! exerts 
at the point r on another ring of mass m: 

*M = -/^ (U 

where G is the gravitation constant, p' is the mass density of the differential volume element 
dV , A is the separation between the perturbing mass element p'dV at r' and the field point 
r, and the integration proceeds over the three dimensional extent of ring m'. In the following, 
primed quantities will refer to the perturbing ring m! and unprimed quantities will refer to 
the perturbed ring m. Each ring can be thought of as a swarm of numerous particles all 
having a common semimajor axis a and an identical mean orbital eccentricity e, inclination 
i, longitude of periapse u, and longitude of ascending node Q. It is also assumed that these 
particles have an isotropic dispersion velocity c that gives rise to the ring's finite radial as 
well as vertical half-thickness h ~ c/n, where n is the ring's mean motion. It shall also 
be assumed that the density p' varies only in the azimuthal direction due to the keplerian 
motion of the ring's particles; in this case the density is p' = X' /Ah' 2 where 



m'r' 



\' = ""' (2) 
2na' 2 Vl^7 2 

is the ring's azimuthal mass per unit length (Murray & Dermott 1999). In cylindrical 
coordinates r'= (£', <fi', z') and dV = £'d£'d(f)'dz f . If the slight radial variations in equation 
(l)'s integrand are ignored (i.e., J i'di' ~ 2r'h'), the potential becomes 

where z' (<f>') is the longitude-dependent height of the perturbing ring's midplane from the 
z = plane. Of course the perturbed ring m also has a radial and vertical half-thickness h, 
and it is useful to form an effective potential by averaging $'( r ) over ^ ne radial and vertical 
extent of ring m: 

^^^L^L 2k^ r)s -L^— Q (4) 

where 

[ z ° +h dz f z '° +h ' dz' r 
Q ^J Zo _ h 2hJ z ,_ h , 2VA (5) 
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where again the slight variations in the integrand with £ are ignored, and that only A is 
assumed to be sensitive to the variations in z and z' . Note that this averaging of $' is 
essential in order for the algorithm developed below to conserve angular momentum. 

The next task is to evaluate the double integrals in Q. First note that the separation A 
between the perturbing mass element p'dV at r' and the field point r obeys A 2 = r 2 + r' 2 — 
2y?(r 2 - z 2 )(r' 2 - z' 2 ) cos(0' - cf>) - 2zz' where r 2 = £ 2 + z 2 . Setting a = r'/r, (3 = z/r, and 
f3' = z'/r'i then 

~l + a 2 -2acos{cf)' -cf)) + {(3 2 + (3' 2 )acos{cf)' - cf) - 2af3f3' (6) 

to second order in the /3's, which are of order the rings' inclinations and are assumed small. 
Inserting this into the expression for Q yields 



Q= dp df3' I + a 2 -2a cos Acf) + (f3 2 + f3')a cos Acf)- 2af3f3' 



7_^/ 4{)()V« cos A0 



(7) 



where f) = /i/r ~ /i/a and ()' = h'/r' ~ /i'/a' are the fractional half-thicknesses of rings m 
and m', /3 G = z D /r and = ^/r' are the rings' midplane latitudes, Acf) = <f>' — <f>, and the 
right-hand side of equation (7) is the result of doing the integration in (3, where 

f3'a — f3 Q a cos Acf> — fya cos Acf) — a/ a cos Acf>(T + e) 

T = = (8) 

f3'a — f3 a cos Acf> + \)a cos Acf> — y a cos Acf>{T — e) 

with r = 1 + a 2 - 2a cos A(f) + ((3 2 + t) 2 + (3' 2 )a cos Acf) - 2a/3'f3 and e = 2fj(/3 cos A0 - /3')a. 
The /3's and the f)'s are assumed small, so £ is second order in the small quantities and is 
negligible when compared to other terms. Thus 



1 — [13' a — f3 Q a cos Acf) — t)a cos Acf)] / v 'a cos AcpT la cos Acf) 

~ 1 - [f3'a - f3 D a cos Acf) + fja cos Acf)]/ ^a cos AcpT ~ + V f ^ 

so In T ~ 2\)\/ a cos Acp/T. This is inserted back into equation (7) and the remaining integral 
over f3' is evaluated similarly, yielding Q ~ ln(A)/2h'\/a; cos Acf) where 

1 — [f3 a a — (3' a cos Acf) — tya cos Acf)] /Jocos A0(^ + Z + £) . la cos Acf) 

A = ^^^^^^=^^^^^^= ~ 1 + 21) ■ 



1 - [f3 a - f3' Q a cos Acf) + fj'a cos A0]/ cos A0(^ + Z - f ) V * + 2 

(10) 



- 8- 



with = 1 + a 2 - 2a cos A0(1 - H 2 ) ~ (1 + a 2 ) (1 + # 2 ) - 2a cos A0, # 2 = (f) 2 + f) ,2 )/2, 
Z = (/3 2 + ^ 2 )o;cos A0 — 2a/3 (3' and £ = 2f)'(/^cos A0 — (3 Q )a is another negligible term. 
Inserting In A ~ 2ty^a cos A0/(\l/ + Z) back into Q yields 



1 ~ #-V2 _ 1^-3/2 

+ 2 2 



;n) 



A Fourier expansion of ^ s will be useful, i.e., \1/ s = | S 



s ~ I Em=-oo ^ cosm(0'-0) where 



cos(m<j))d<f) 



7cJ {(l + a 2 )[l + i(l) 2 + r) - 2a cos 0} 



(12) 



is the softened Laplace coefficient. The usual unsoftened form is 6^ (a) = b~i m \a, 0, 0). 
These two coefficients are nearly equal when a is far from unity, but the softened form is 
finite at a — 1 whereas the unsoftened form diverges. Writing Q in terms of softened Laplace 
coefficients thus gives 



oo 



9 ^ cosm(0' - 0) <j //,":,' 

m=— oo 



^(/5 2 + ^ 2 )acosA0-^> 



J 3/2 



which is then inserted in equation (4) to get the perturbing ring's potential 



< $ 



(/3 2 + / ^)acos(0'-0)-a/y^ 



(13) 



°3/2 

(14) 



The final task of this Section is to write the rings' coordinates in terms of its orbit 
elements 



r ~a(l — e cos z/ e 2 H — e 2 cos 2v) 

=u + ~ a) + M + 2e sin M 



sini sin(0 — fi) 



(15a) 
(15b) 
(15c) 



where v is the true anomaly of a ring-element at r and M = nt is the corresponding mean 
anomaly. Inserting equations (15) into the potential < $' >, expanding to second order in 
the e's and i's, doing the 0' integration in equation (14), and time-averaging the resulting 
expression over the orbital period of ring m then yields the time-averaged potential < $' > 
experienced by ring m due to ring m' assuming small eccentricities and inclinations: 



Gm! 



l -b^ 2 + l -{e 2 + e' 2 )f + \ee'co^'-^)g 



+ f)ab$ 2 + i»' cos(fi' - n)a6« 



(16) 
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where the / and g functions are 

f(a, h, h') = [2a^ + a 2 ^j 6® = ab^ 2 - 3a 2 ff 2 (2 + H 2 )bf} 2 (17a) 
(/(a, h, h') = (2 - 2a A - a 2 6« = -a&$ + 3a 2 iJ 2 (2 + tf 2 )^ (17b) 

where if 2 = |(f) 2 + h /2 ) and a has been redefined as a = a' /a. The right-hand side of 
equations (17) are derived in Appendix A, and it is shown in Appendix B that the softened 
Laplace coefficients can be written in terms of complete elliptic integrals. Consequently, 
functions /, g, and all can be rapidly evaluated without relying upon a numerical 
integration of equation (12). It is also noted that the b~i m \ f, and g functions obey the 
following reciprocal relations 

b^\a-\^^)=aH^ ] M^) (18a) 

/(a-\b',fO=a/(a,f),f)') (18b) 
g(a- 1 ,t)',t))=ag(a,l),l)'). (18c) 

These relations will be used in Section 2.2.1 to show that the equations of motion developed 
below conserve angular momentum. 

The laborious procedure of expanding, integrating, and then time-averaging < $' > 
is not included here since a similar analysis can be found in Murray & Dermott (1999) 2 . 
Since the terms proportional to e' 2 and i' 2 as well as the first term in equation (16) do not 
contribute to the resulting dynamical equations, they may be neglected. The disturbing 
function R for ring m due to another ring w! is — lx the surviving terms in < $' >, i.e., 



-/e 2 + -^gee cos(c/ — uj) — -ctb^ 2 i 2 + -ab^j 2 ii' cos(f2' — Q) 



(19) 



Note that when h = h' = the disturbing function for a point mass m perturbed by point 
mass m! is recovered (Brouwer & Clemence 1961; Murray & Dermott 1999), which is to be 
expected since both a point mass as well as a thin ring have the same disturbing function to 
this degree of accuracy (Murray & Dermott 1999). 

It should also be noted that many celestial mechanics texts develop two distinct expres- 
sions for the disturbing function R, one due to a perturber in an interior orbit having a' < a, 



2 Actually Murray & Dermott (1999) derive the time-averaged acceleration (rather than a potential) that 
ring ml exerts on to, which they insert into the Gauss equations to obtain a set of dynamical equations 
equivalent to that obtained here when f) = 0. 
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and another expression due to an exterior perturber having a' > a [c.f., Brouwer & Clemence 
(1961); Murray & Dermott (1999)]. However this pairwise development is unnecessary in 
this application since equation (19) is valid for a = a' /a < 1 as well as for a > 1. Indeed, it is 
straightforward to show that these pairs of disturbing functions, such as equations (7.6) and 
(7.7) in Murray & Dermott (1999), are in fact equivalent to equation (19) with \) = = 0; 
they only appear distinct since one is a function of a and the other is actually a function of 



a 



In terms of the variables 



h —e sin Q 
k —e coscu 



p =i sin Q 
q —% cos Q, 



(20a) 
(20b) 



the disturbing function can then be written 



R = n a 



2„2 



m 



M Q + m 



l f{h 2 + k 2 ) + -g{hh' + kk') - \ab§ 2 (p 2 + q 2 ) + 



8 



8 

\<^i] 2 {pp' + qq) 



(21) 



where the mean motion n = sjG(M & +m)/a? with M being the solar mass. 



2.1.1. the N ring problem 

For the more general problem of N perturbing rings, replace the perturbing ring mass 
m' with mk and give all the other primed quantities the subscript k. The disturbing function 
R — > Rj for the perturbed ring having a mass m — > rrij is sum of equation (21) over all the 

other rings k ^ j. Setting <x, fc = a^ja^ rij = G(M Q + rrij) / a|, and 

Ajk = f ( mJ+wJ ^ ^ w j ^ fc (22b) 

^' = S ( M ^ m . ) a 3kb§ 2 (a jk , f) fc ) (22c) 
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that sum becomes 



-Ajj{h) + fcj) + A jk(hjh k + kjk k ) + -Bjjipj + gj) + Y B jk (p jPk + qj q k ) 

k+j ~ k^j 

(23) 

in the notation of Murray & Dermott (1999). The A jk and the B jk can be regarded as two 
N x N matrices A and B whose entries describe the magnitude of the mutual gravitational 
interactions that are exerted among the N rings. In the following, quantities having a j 
subscript will always refer to the perturbed ring in question while the k subscript will always 
refer to another perturbing ring. 

The time-variation of the rings' orbit elements is given by the Lagrange planetary 
equations; to lowest order in the e's and i's these are (Brouwer & Clemence 1961; Murray & 
Dermott 1999) 

dhj dRj/dkj 



J2 A 3kh (24a) 



dt rijO 2 - 

3 3 k=l 

dkj dRJdhj . , . 

J J k=l 

(24c) 

3 3 k=l 



dqj dRj/dpj V- R „ (9Ar]] 

3 j k=i 

and their well-known Laplace-Lagrange solution is 

N 

hj (t) = E ji M^t + A) (25a) 
i=i 

N 

kj(t) = Y E ji c os{g t t + Pi) (25b) 
i=\ 

N 

Pj (t) =J]/ ji sin(/ i * + 7 i ) (25c) 
%=\ 

N 

% (*) = Y ^ cos (^ + 7j ) ( 25d ) 

1=1 

where gi is the i th eigenvalue of the A matrix, fi is the i th eigenvalue of B, is the N x N 
matrix formed from the N eigenvectors to A and Iji is the matrix of eigenvectors to B, and 
(3i and 7j are integration constants. 
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To apply this rings model, first assign to the N rings their masses rrij, their semimajor 
axes a,, and their fractional half-thickness Planets are represented as thin rings having 
f)j = 0. Then construct the system's A and B matrices and compute their eigenvalues Qi, 
fi and the eigenvector arrays and Iji. The rings' initial orbits ej(0), ij(0), %(0), and 
Qj(0) are then used to determine the integration constants and 7« as well as to rescale 
the eigenvectors Eji and Iji such that equations (25) agree with the initial conditions. A 
handy recipe for this particular task is given in Murray & Dermott (1999). Equations (25) 
are then used to compute the system's time-history, and orbit elements are recovered via 
e 2 = h 2 + i 2 = p 2 + q 2 , tancj, = hj/kj, and tanfij = Pj/qj. 

The rather laborious derivation given above thus confirms the assertion by Tremaine 
(2001) that one only needs to soften the Laplace coefficients in order to use the Laplace- 
Lagrange solution for calculating the secular evolution of a continuous disk. However Section 
2.2.1 will show that this softening must be done judiciously such that equation (18a) is obeyed 
in order for the solution to preserve the system's angular momentum. 

2.2. tests of the rings model 

Several tests have been devised in order to demonstrate that the rings model described 
above behaves as expected. 

2.2.1. angular momentum conservation 

The equations developed above conserve the system's total z-component of angular 
momentum L z = mjnja 2 ^jl — e 2 cos(ij). To show this, expand L z to second order in the 
e's and i's, which is the same degree of precision to which the disk potential is developed 
above. This gives L z ~ L — L e — Li where 

L e = l - "'.,»j"j<j ( 26a ) 

j 

'■' i Yl '".i"j"Vj ( 26b ) 

j 

and L = mjUja 2 is a constant. It is shown in Appendix C that the dynamical equations 
(24) conserve angular momenta to second order in the e's and i's, that is, dL e /dt = and 
dLi/dt = 0. When the rings model is used to calculate the secular evolution of a 'sparse' 
system, such as the Jupiter-Saturn system described below or one composed of all four giant 
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planets, the angular momenta L e and Lj are preserved to near machine limits to a fractional 
precision of ~ 10~ 7 in this single floating-point implementation. This is true regardless of 
whether the planets are represented as thin rings having fyj = or as thick rings with > 0. 
However L z conservation is a little worse, ~ 10~ 6 , in the 'crowded' systems described in 
Section 4 that consist of a few planets plus a disk comprised of many closely-packed rings. 
However these errors are always smaller, by a factor of ~ 10 4 to 10 5 , than the angular 
momenta associated with the disturbances and waves seen in these disks. Thus the wave- 
like disk behavior reported below is real and is not due to a diffusion of some numerical 
error. 



2.2.2. Jupiter and Saturn 

Murray & Dermott (1999) use the Laplace-Lagrange planetary solution, equations (25), 
to study a system composed of Jupiter and Saturn. The rings model developed here repro- 
duces that system's A and B matrices, eigenvector arrays and eigenvalues gi and 
fi, and integration constants $ and 7«, to the same precision quoted by Murray & Dermott 
(1999), provided one sets nj = GM Q / and \)j = 0. The rings model also reproduces 
the figures in Murray & Dermott (1999) that show this system's orbital history, as well as 
the figures showing the forced orbit elements of numerous other massless 'test-rings' orbit- 
ing throughout this system. However is should be noted that rigorous angular momentum 
conservation actually requires setting nj = JG{M Q + rrij)/a^ instead. 



2.2.3. precession in an axisymmetric disk 

Heppenheimer (1980) points out that a massless test particle orbiting in a smooth, ax- 
isymmetric disk experiences a regression of its longitude of periapse, i.e., Cj < 0, which is the 
opposite of the usual prograde apse precession that occurs throughout the Solar System. As 
Ward (1981) shows, it is the nearby disk parcels whose orbits actually cross the test particle's 
orbit that drives periapse regression at a rate that exceeds the prograde contribution from 
the more distant parts of the disk. 

The rings code also reproduces this phenomenon. An annulus having a surface density 
a', mass 5m' = 2-ira'a'da', a semimajor axis a', and a fractional half-thickness h' contributes 
5R = n 2 a 2 [f(a, 0, h')e 2 — 0:63*2(0:, 0, t)')i 2 ]5m' /8M to the particle's disturbing function [see 
equation (19) with e' — i' — 0]. Adopting a power-law in the disk surface density, a' = 
a(a)a~ r where a = a' /a, the total disturbing function integrated across a semi-infinite disk 
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j 5R= ^ d n 2 a 2 [he 2 - I n t 2 ] (27) 
where fid = ncr(a)a 2 /M Q = nGa/an 2 is the so-called normalized disk mass and 
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h =^ / a 1-r /(a, 0, (28a) 



o 



1 r°° 

/n=-y a 2 - r b^ 2 (a,0,t)')da. (28b) 

Note that I& and Iq are double integrals according to the definitions of and /, equations 
(12) and (17a). However these integrals are analytic for selected power laws r. For instance, 
setting r = 1 or r = 2 and then instructing symbolic math software such as MAPLE to do the 



radial integration first and the angular integration second yields I& = —1/ yl + f) /2 /2 ~ — 1 
and I Q = [f)' _1 + ^1 + f) ,2 /4 + fj'/^/v^ + (j'74)(l + f) /2 /2) ~ 1/fj' for small fj'. The 
Lagrange planetary equations then give the test-ring's precession rates: 

- dR / 9e r 

uj ~ — = Icj(j, d n ~ -/i d n (29a) 

na z e 

« = ~ -/W)'- (29b) 

nan 

Very similar precession rates were previously derived in Heppenheimer (1980) and Ward 
(1981). Note that the model rings start to overlap when their fractional radial thickness 2(j' 
exceeds the rings' fractional separation 5 = Aa/a, so a massless test ring should precess at 
the above rates when the disk-rings are sufficiently overlapping. 

These expectations are tested by constructing a 50 M e disk having an r = 2 power-law 
surface density using 200 circular, coplanar rings arranged over 10 < a < 100 AU, with a 
number of thin, massless 'test-rings' also orbiting within this disk. Several simulations are 
performed with the massive rings having a variety of thicknesses ()'. As expected, those test- 
rings that reside far from the disk's edges precess at rates given by equation (29) whenever 
the disk-rings are sufficiently overlapping, namely, when [)' > 25. 



2.2.4- precession in an eccentric disk 

Although a particle's periapse will experience retrograde precession when embedded 
in an axisymmetric disk, prograde precession is possible in a non-axisymmetric disk. For 
instance, prograde precession is evident in the N-body simulations of an eccentric stellar disk 
orbiting the putative black hole at the center of the Andromeda galaxy M31 (Jacobs and 
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Sellwood 2001). These simulated disks have masses ~ 0.1 times the central mass, with the 
interior parts of the disk being progressively more eccentric. These simulations reveal a long- 
lived overdense region in the direction of apoapse that persists due to a coherent alignment 
of the particles' periapse, with the density pattern rotating in a prograde sense at rates that 
increase with the disk mass. It should be noted that the disk particles' eccentricities are 
not always small, so these N-body simulations cannot be used as a quantitative benchmark 
for the rings code. Nonetheless, it is comforting to find that the rings code does indeed 
reproduce the density patterns seen in the N-body disks that precess at rates very similar 
to that reported in Jacobs and Sellwood (2001). 



Section 4 will use the preceding rings model to demonstrate that apsidal density waves 
and nodal bending waves may have once propagated throughout the Kuiper Belt. An apsi- 
dal wave is a one-armed spiral density wave that slowly rotates over a periapse precession 
timescale. Similarly, a nodal wave is a one-armed spiral bending wave that rotates over a 
nodal precession timescale. 

A brief review of spiral wave theory is in order. Many of the waves' properties, such 
as their wavelength and propagation speed, are readily extracted from the waves' dispersion 
relations. These dispersion relations are usually obtained from solutions to the Poisson and 
Euler equations for the disk. However the following will show that these dispersion relations 
may also be derived from the Lagrange planetary equations. 



The disturbing function R(a) for the disk material orbiting at a semimajor axis a is 
obtained from equation (19) with the perturbing mass m! replaced by the differential mass 
dm' whose contributions are integrated across a semi-infinite disk: 



3. Spiral Wave Theory 



3.1. Apsidal Density Waves 




a 



l-r 




where e'(a) and Co 1 {a) are the orbit elements of the perturbing parts of the disk, the unprimed 
quantities refer to the perturbed annulus at a, r is the power-law for the disk's surface 
density variation, and a constant fractional thickness f) is assumed throughout the disk. The 
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Lagrange planetary equations then give the disk's periapse precession rate at a: 

■ . . dR/de j . . . 

00(a) ~ =— = n d n{Iu + hw) (31) 

na z e 

where — — 1 is from the left term in equation (30) and is the contribution by an undis- 
turbed disk to its own precession [see equations (28-29)], and the right term is 

1 f°° e' 

Idw = - - cos(^' - od)^ r g(a, \), i))da, (32) 

2 Jo e 

which is the relative rate at which the density wave drives its own precession. 

It is expected that the eccentricities associated with a spiral density wave will vary only 
slowly with distance a such that e'(a)/e ~ 1. The spiral wave will also organize the disk's 
longitude of periapse Co such that it varies as oo(a,t) = Cot — f a k(A)dA, where k(a) is the 
wavenumber and A = 2n/\k\ is the radial wavelength. If the spiral wave pattern is tightly 
wound such that A < a and \ka\ ^> 1, then the dominant contributions to I dw are largely 
due to the nearby parts of the disk where a = a' /a ~ 1 ± A/a and Co' — Co = — k(A)dA ~ 
—k(a' — a), while the contributions from the more distant parts of the disk tend to cancel 
owing to the rapid oscillation of the cosine factor. Thus we can set a = 1 in equation (32) 
except where it appears as the combination x = a — 1 where \x\ <^ 1. In this case the 
softened Laplace coefficients that are present in the g function can be replaced with the 
approximate forms that are valid for \x\ <^ 1 [see equation (B5)], so g becomes 

2 2[) 2 -x 2 

9{X) " n WTxT ' (33) 

It is also permitted to extend the lower integration limit in equation (32) to — oo in the 
tight-winding limit, so 

I dw ~- cos(\ka\x)g(x)dx = \ka\e~^ ka \. (34) 

2 J -OO 

And finally, if the wave is to remain coherent across this disk, this self-precession must 
occur at the same constant rate to throughout the disk, so Cb(a) = to ~ fx d (\ka\e~^^ ka ^ — l)n. 
Note that this disturbing frequency oo, which is also called the pattern speed, can also be 
identified as any one of the eigenfrequencies gi that appear in equations (25). Usually it 
is another perturber that is responsible for launching the wave and causing the disk to 
precess in concert at the rate oo, and this is usually at a rate that dominates over the non- 
wave contribution to the disk's precession, i.e., \oo\ ^> jj, d n. This then yields the dispersion 
relation for tightly-wound apsidal density waves: 

oo ~ e-^ ka ^ d \ka\n (35) 
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which has the dimensionless form 



u e (K) = Ke 



-K 



(36) 



where u e = \[2\)u j \i d n is the dimensionless frequency and K = \/2t)\ka\ the dimensionless 
wavelength; Figure 1 plots ui e versus K. A more general dispersion relation for an m-armed 
spiral wave in a stellar disk is given in Toomre (1969), and in the limit that \uo\ <C n the 
resulting formula with m — 1 behaves qualitatively quite similar 3 to equation (36). Note 
that equation (35) also recovers the usual dispersion relation uj = /i d \ka\n for apsidal waves 
in an infinitesimally thin disk when h = 0. 

Note that ui e (K) > and that it also has a maximum at K — 1 where uo e (l) = exp(— 1). 
Since oo e is a function of semimajor axis a, the restriction < uo e (a) < 0.368 indicates that 
apsidal waves can only propagate in a restricted interval in a. This restriction can also be 
viewed as a constraint on the disk thickness, namely, 



Alternatively, equation (37) can also be viewed as an upper limit on the frequency u or a 
lower limit on the disk mass /id wherein wave-action is permitted. 

Waves having a wavenumber K < 1 are called long waves since they have a wavenum- 
ber \ki\ — (cu /n) / fj,da and the longer wavelength Al = 2ir/\k,L,\ — 2ir/id(n/ui)a, while 
short waves have wavenumbers K > 1 or \ks\ > l/y/2t)a and the shorter wavelength 
\ s = 2ir/\ks\ < 2v / 27r()a. These long and short density waves also correspond to the g 
and p modes, respectively, of Tremaine (2001). Writing Al in term of [)q and then requiring 
[) < f)g also means that apsidal wave can propagate wherever A^ > 24h. 

The rate at which apsidal waves propagate across the disk is given by their group velocity 



where Sk = sgn{k). Waves with s fc = +1 are called trailing waves, and their longitude of 
perihelia Cj decreases with increasing semimajor axis a, while Sk = —1 are leading waves 
whose longitude increases with a. Since the waves' group velocity c g is proportional to the 
slope of the w e (K) curve, the site where K — 1 is a turning point where wave reflection 
occurs; in galactic dynamics this reflection site is known as a Q barrier. A long wave that 
approaches the Q barrier from the K < 1 side of Figure 1 thus reflects as a short wave as it 



3 In this limit, Toomre's dispersion relation becomes lot{K) = KT where the reduction factor T is a more 
complicated function of K. However a numerical evaluation of this function shows that Wt{K) has the same 
form as the ui e (K) curve shown in Figure 1. 



f) < 0.260^ d |n/w 



(37) 




(38) 



continues along the K > 1 side of the curve. The simulations shown in Section 4 also show 
that long trailing waves that instead strike a disk edge also reflect as short trailing waves. 



3.2. surface density variations 

The compression or rarefaction among the disk's rings, or streamlines, is (dr/da)^ 1 , 
which is also be the relative change in the disk's surface density a ja associated with a density 
wave [c.f., Borderies, Goldreich & Tremaine (1985)]. For a small amount of compression, 
a /o- = 1 + Aa/a Q where the fractional surface density variation is 

Aa /9r\ _1 d(ea) ,, _ N dtu . . . /nnS 

-=(*) -l = -^o M («- W ) + « s «.(#-«) (39) 

to lowest order in e. The second term dominates over the first in the tight-winding limit so 
\Aa ja \ ~ 0(e\ka\). Density waves are nonlinear when \Aa/a a \ > 1, and these large density 
variations are a consequence of overlapping streamlines. For long density waves having a 
wavenumber \ki\ = u/fiddn, the disk's streamlines will cross when the waves' eccentricities 
exceed ~ fi d u/n, while streamline crossing occurs among short waves when eccentricities 
exceed es ~ V^f). As the simulations of Section 4 will show, a dynamically cool Kuiper 
Belt is very susceptible to the propagation of short nonlinear density waves that facilitate 
streamline crossing. Depending upon the relative velocities of these crossed streamlines, 
apsidal wave-action might encourage accretion or else enhance collisional erosion among 
KBOs. 



3.3. Nodal Bending Waves 



The derivation of the dispersion relation for nodal bending waves, and its analysis, 
proceeds similarly. The disk's integrated disturbing function is 



R(a) = --fi d n 2 a 2 J a 2 ~ r bf} 2 {a, f), Jj) 
so the Lagrange planetary equation gives 



h 2 - %i! cos(Q' - O) 



da. 



Q(a) ~ dR l® 1 = jj, d n(I bw - I n ) 



where Q{a, t)=(lt- f k(A)dA and 

Ibw 



a 2 ~ r b W 



3/2 ( a ) hi h) cos[A;a(a — l)]da ~ — -^=- — 



(40) 



(41) 



(42) 
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is the bending wave's contribution to its own precession. The contribution from the undis- 
turbed disk is Iq ~ l/-\/2f), with the additional \[2 factor the result of changing the middle 
argument in equation (28b) from to f). Since u = Cl is a constant for a coherent wave, the 
dispersion relation for tightly-wound nodal bending waves is 

u ~ _g_ ( e -^i*«i _ n _ (43) 

Note that the usual dispersion relation for nodal waves in an infinitesimally thin disk, u ~ 
— /id\ka\n, is obtained when rj — > 0. 



The dimensionless form of the dispersion relation is 

Ui(K) = e- K - 1, (44) 

where = y/2t)u / fian, which is plotted in Figure 1. As the Figure shows, nodal bending 
waves can propagate only in regions where — 1 < Ui(a) < which similarly limits the disk 
thickness to f) < 2.72f)Q. The nodal waves' group velocity is 

dco 

which indicates that nodal waves tend to stall, i.e., |c 5 | — > as they approach the Ui — — 1 
boundary. Note that this dispersion relation only admits a long wavelength solution having 
a wavenumber \ki\ ~ —u/fidna and a wavelength Al ~ 2njid\n/oj\a for waves far from the 
stall-zone. Since [) < 2.72F)q, the disk can sustain nodal waves wherever > 9h. But if an 
outward-traveling long leading wave with s k = — 1 encounters a disk edge, it will reflect as 
an Sk = +1 long trailing wave. Examples of nodal wave reflection and stalling are also given 
below. 



4. The Secular Evolution of the Primordial Kuiper Belt 

Using the recipe given in Section 2.1.1, the rings model is used to compute the secular 
evolution of the primordial Kuiper Belt as it is perturbed by the four giant planets. In these 
simulations the giant planets are represented by thin [) = rings whose initial orbits are 
their current orbits, while the Kuiper Belt is represented by 500 rings whose semimajor axes 
extend from 36 AU out to 50-70 AU. The location of the Belt's inner edge is chosen such 
that only the outermost radial and vertical secular resonances, the u 8 and the ui 8 , reside in 
this disk near 40 AU when of low mass. The semimajor axes of each Belt-ring increases 
as Oj+i = (1 + 8)a,j where the rings' fractional separation 5 is typically ~ 0.001. The rings 
fractional half-thickness I) is always in excess of 25, as is required to get the correct apse 
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precession in an axisymmetric disk (see Section 2.2.3). The Belt-rings initial eccentricities 
and inclinations are zero, with all inclinations being measured with respect to the system's 
invariable plane. The mass of each ring is chosen such that the Belt's surface density a(a) 
varies as a~ L5 . For this configuration, if the total Kuiper Belt mass over the 36 < a < 70 AU 
zone is M to tai, then the total mass in the 'observable' 30 < a < 50 AU zone that is presently 
accessible to astronomers would be Mkb = 0.67M tota i had the above surface density law 
extended inwards to 30 AU. For these systems the normalized disk mass is fi d ~ M KB /M Q . 



4.1. A M KB = 10 M ffi Example 

Figure 2 shows a snapshot of apsidal density waves as they propagate across an b — 
10 M® Kuiper Belt having a half-thickness f) = 55 = 0.0067. Since f) ~ 0.2\jq, the necessary 
disk-conditions for the propagation of apsidal and nodal waves are well-satisfied. Initially, 
a long trailing density wave is launched at the Belt's inner edge. This wave is really more 
like a pulse ~ 5 AU wide, and Figure 2 shows that by time t = 2 x 10 6 years the wave has 
just started to reflect at the disk's outer edge at 70 AU. The greyscale map of the disk's 
surface density variation, Aa/a, is obtained using equation (39), and this map also provides 
a historical record of this system's wave-action. It should be noted that equation (39) is 
quantitatively correct only when \Aa/a\ <C 1, a condition that is rarely satisfied by the 
results obtained here. Nonetheless, equation (39) is still useful in a qualitative sense since it 
will indicate the disk-regions where large surface density variations as well as orbit-crossing 
can be expected. As the outer edge of the Aa/a map shows, the main apsidal density 
wave-pulse at 67 < a < 70 AU has just reflected at a short trailing wave, and this nonlinear 
wave having \Aa/a\ > 1 will completely dominate the disk's appearance at later times 
as it propagates inwards. But until this happens, the bulk of the disk's appearance over 
45 < a < 67 AU is still dominated by lower-amplitude long waves that are following behind 
the main density pulse. Note also that short leading waves were also emitted at the disk's 
inner edge, but as the density greyscale shows, they have only propagated out to a = 45 AU 
thus far owing to their slower group velocity [see equation (38)]. These short waves are well- 
resolved in the sense that their radial wavelength of As ~ 1 AU spans about 15 disk-rings. 
Although these short waves are seen in the e(a) and Cj{a) plots as only tiny wiggles over 
36 < a < 45 AU that are superimposed on top of the disk's longer-wavelength behavior, 
the density greyscale shows that the short waves dominate the inner disk's appearance. The 
online edition of Fig 2 is also linked to an animated sequence of these figures that give the 
system's complete time-history. That animation shows that by time t ~ 1 x 10 7 years, 
nonlinear short waves will have swept across the entire disk, and they result in large surface 
density variations \ Aa/a\ > 1 over radial wavelengths of As ~ 1 AU. 
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Figure 2 also shows that by time t = 2 x 10 6 years, a long leading nodal bending 
wave-pulse has already propagated across the disk where it has reflected at the disk's outer 
edge and just started to return inwards as a long trailing bending wave. But when that 
pulse reaches the disk's inner edge, a portion of the wave's angular momentum content will 
continue to propagate further inwards where it will give a small kick to the giant planets' 
orbit elements, while the remaining wave-pulse reflects again and propagates outwards. The 
same phenomenon also occurs among the apsidal density waves. Thus after a few reflections, 
a single wave-pulse will lose its initial spatial coherence by spawning multiple wave-trains 
that, in this friction-free model, forever roam about the Belt. This ultimately results in a 
rather wobbly-looking standing density wave pattern that varies over the short wavelength 
scale \s ~ 1 AU, as well as a standing bending wave pattern that varies over a much longer 
scale \ L ~ 40 AU. 

A dynamical 'spectrum' of this Kuiper Belt is shown in Figure 3, which plots all of the 
eccentricity eigenvector elements \Eji\ for each of the disk-rings versus their eigenfrequency 
gi, as well as inclination eigenvector elements \Iji\ versus eigenfrequency The upper 
figure is quite reminiscent of the findings of Tremaine (2001), who showed that a gravitating 
disk tends to exhibit its strongest response to slow radial perturbations via modes having 
discrete patterns speeds uj that can be identified with any of the peak frequencies seen in 
Figure 3. Figures such as these are quite useful for identifying the patterns speeds uj that 
are associated with the density and bending waves that propagate in the disk. 



4.2. Variations with Kuiper Belt Mass M KB 

Simulations have been performed with Kuiper Belts having masses M^b — 0, 0.08, 
0.2, 2, 10, and 30 M e in the observable 30 < a < 50 AU zone (again this would be the 
Belt's mass assuming its surface density a (a) oc a~ 15 were to extend solely over that region). 
The results are summarized in Figure 4 which shows the maximum eccentricities e max and 
maximum inclinations i max achieved by the rings in each simulation. The Belt's radial width 
is indicated by the breadth of the curves in Figure 4, which ranges from 34 AU in the 
higher-mass Belts to 14 AU for the Mkb = 0.08 M e system. Each simulation uses 500 
disk-rings having a fractional half-thickness f) = 25, so the three higher-mass systems have 
f) = 0.0027 while the lower-mass disk M^b — 0.2 M e is somewhat thinner with {) = 0.0015 
and the Mkb — 0.08 M e system has f) = 0.0010. The rings in these lower-mass disks 
are more closely packed so that their shorter-wavelength density waves are well resolved, 
and they are also made thinner so as to push their Q-barriers a bit further downstream. 
The disk's initial orbits are e = = i, except for the M KB = system which instead 
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adopts the forced orbit elements appropriate for a massless disk [e.g., Brouwer & Clemence 
(1961); Murray & Dermott (1999)]. These systems are evolved until their initial density and 
bending wave pulses have had the opportunity to reflect multiple times and have largely 
dissolved into standing wave patterns. The lower-mass disks necessarily have longer run- 
times due to their waves' slower propagation speeds [see equation (38)]; these run-times are 
1 x 10 9 years, 2 x 10 9 years, 1 x 10 8 years, 2 x 10 7 years, and 1 x 10 7 years, respectively, 
for the M KB = 0.08, 0.2, 2, 10, and 30 M e systems. Once the standing wave pattern has 
emerged, the rings' instantaneous eccentricities and inclinations range over < e < e max and 

< l < W- 

As Figure 4 shows, there are no peaks in the e max and i max curves for the higher mass 
disks having M KB > 2 M e , indicating that there are no secular resonances in the disk itself; 
any such resonances likely lie between the orbits of Neptune and the disk's inner edge at 36 
AU. The simulations of these higher-mass disks show long trailing density waves and long 
leading bending waves being launched at the disk's inner edge. These waves sweep across 
the disk, reflect at the disk's outer edge, and return as short trailing density waves and long 
trailing bending waves, similar to the history described in Section 4.1 and seen in Figure 2. 

However secular resonances at a ~ 41 AU are quite prominent in the lower-mass disks 
having Mkb = 0.08 and 0.2 M®. These are sites that launch long density and bending waves, 
and as Figure 2 shows, the density waves are able to propagate out as far as a = 45 and 49 
AU, respectively, where they are reflected by a Q-barrier and return as short density waves. 
These reflection sites occur where f) ~ 2.2f)g and () ~ l-4hg, which shows that equation 
(37) is an approximate yet useful indicator of where apsidal wave are allowed to propagate. 
Note also the large amplitudes of the density waves in these low-mass disks, 0.3 < e max < 1, 
which clearly violates the model's assumption of low eccentricities. Thus these particular 
curves should not be taken at face value. Nonetheless, they do indicate that apsidal waves 
in a low-mass Kuiper Belt may result in large eccentricities, and Figure 2 shows that sizable 
inclinations can also result due to nodal bending waves. Indeed, in the Mkb = 0.08 M e 
disk, maximum inclinations are typically i ~ 20°, which is comparable to the Main Belt's 
high-inclination component. 

The larger e's and i's seen in the lower-mass disk's are a consequence of the giant 
planets transmitting a small fraction of their initial angular momentum deficit 4 (AMD) into 
the disk in each simulation. In each simulation, the giant planets deposit ~ 0.005L e and 
~ O.lLj into the disk's inner edge, which waves then transport and smear out across a vast 
swath of the Kuiper Belt. Since the AMD deposited in the disk is roughly constant in 



e.g., the L e and Li of equations (26). 
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each simulation, the lower-mass disks exhibit larger e and % excitations. It should also be 
noted that computational limitations in dynamical studies of the Kuiper Belt, particularly 
N-body simulations, are often require treating the Belt as a swarm of massless test particles. 
However the comparison of the M^b — curve to the M^b > curves in Figure 2 shows 
that the endstate of a Kuiper Belt having even just a modest amount of mass can be radically 
different from one naively treated as massless. 

4.3. Variations with Disk Thickness f) 

Figure 5 shows how the response of an M^b — 10 M e Belt varies with increasing disk 
thickness f) = (2, 20, 30, 60, 100) x 5 =(0.0027, 0.027, 0.04, 0.08, 0.13). According to equation 
(37), f)g(a) oc a 1 / 2 ~ r= ~ 1 in these r = 3/2 disks, so the Q-barrier will move inwards as the 
disk's f) is increased, as is evident in Figure 5. All of these disks have a normalized disk 
mass fid ~ 3 x 10~ 5 , a mean motion n ~ 0.02 radians/year, and a pattern speed that is 
typically uo ~ 3 x 10~ 6 radians/year, so wave-action is shut off when the disk thickness t) 
exceeds \)q = 0.2Qfi d \n/ 'u\ ~ 0.05. In the thinnest disk with f) = 0.0027, the Q-barrier lies 
beyond the disk's outer edge at 70 AU, so long and then short apsidal density waves are able 
to slosh about the disk's full extent. However the Q-barrier lies in the disk at a ~ 60 AU 
when [) = 0.027 (red curve), at a ~ 53 AU when f) = 0.04 (green curve), and apsidal waves 
are prohibited in the disks with () > 0.08 (blue and purple curves). 

Figure 5 also shows that the nodal bending waves are shut off when the disk thickness 
f) exceeds the somewhat more relaxed criterion 2.72[)q ~ 0.14. Note also the peak at a ~ 53 
AU in the f) = 0.027 disk (red curve) and at a ~ 39 AU in the \) = 0.04 disk (green curve). 
These particular disks admit two spiral patterns, a higher-amplitude spiral that corotates 
with Neptune's dominant eigenmode at the pattern speed of uj ~ — 3 x 10~ 6 radians/year, 
and a lower-amplitude mode that corotates with Uranus at the faster rate uj ~ —1.5 x 10~ 5 
radians/year. Since \)q oc the faster spiral pattern has 2.72f)g ~ 0.03, which is why 

this wave stalls at a = 53(39) AU in the t) = 0.027(0.04) disks while the slower spiral mode 
is able to propagate the full breadth of these disks. 

The behavior of a lower-mass disk with M KB = 0.2 M ffi is shown in Figure 6 for disks 
having fj = (2,5, 10, 15) x 5 =(0.0015, 0.0037, 0.0074, 0.011). A pair of secular resonances 
lie near a ~ 40 AU, and they launch apsidal and nodal waves having pattern speeds u ~ 
±3 x 10~ 6 radians/year. These disks have fid ~ 6 x 10~ 7 and t)Q ~ 0.001, and as the upper 
plot shows, further increases in \) brings the Q-barrier inwards until the wave-propagation 
zone has shrunk down to zero. The lower plot also shows that nodal waves forever slosh 
about the model Kuiper Belt in the f) = 0.0015 disk (orange curve), whereas the peaks in 
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the other curves show that these waves stall at sites ever closer to resonance in progressively 
thicker disks. These Figures also show that when f) 3> \)q, the motions of a non-gravitating 
disk are recovered, namely, the disk's e's and i's are twice the forced motions seen in the 
Mkb = disk (black curves), with the factor of two being a consequence of these disks' 
initial conditions e = = i. 



The primordial Kuiper Belt likely had an initial mass of Mkb ~ 30 M® (see Section 1), 
and accretion models show that the initial KBO swarm must have had dispersion velocities 
less than ~ 0.001 x Keplerian (Kenyon & Luu 1999), so f) < 0.001, [i d ~ 1 x 10~ 4 , and thus 
\)q ~ 0.2 assuming the spiral waves have pattern speeds of \u\ ~ 3 x 10~ 6 radians/year. 
Since () < f)Q, the primordial Kuiper Belt readily sustained apsidal and nodal waves. Figure 
4 shows that in this high-mass environment, these will be rather low-amplitude waves having 
e ~ sini ~ 0.01. These waves will quickly propagate across a Kuiper Belt of width Aa in 
time T prop ~ Aa/\c g \ ~ Aa/^an, so the wave propagation time is 



In the friction-free model employed here, the outward-bound long density waves eventually 
reflect, either at a Q-barrier in the disk (which might lie downstream where f) = \)q) or else 
at the disk's outer edge. The reflected waves then propagate inwards as short density waves, 
and such waves are nonlinear in the sense that their surface density variations Aa /a typically 
exceed unity. Figure 2 shows a snapshot of long and short density waves in a Mkb = 10 M® 
disk. Note that the long waves completely dominate the disk's orbit elements e(a) and uj(a) 
that vary over a wavelength of ~ 10 AU, while the short waves are the tiny variations in 
e(a) and Co(a) that occur over a \ s ~ 1 AU scale at the disk's inner and outer edge. Even 
though the short waves are almost imperceptible in the orbit-elements plots, the greyscale 
map shows that these nonlinear short waves completely dominate the disk's surface density 
structure. 

Figure 4 also shows that the waves in lower-mass systems have higher amplitudes. This 
suggests that wave action may tend to drive disk-planet systems towards an equipartition 
of angular momentum deficits, since the angular momentum content of the waves seen in all 
simulations are ~ 0.5% and ~ 10% of the planets' L e and Lj. 

The large-amplitude waves seen in the M K b = 0.2 M® disk also suggests that apsidal 
and nodal wave action might account for much of the Kuiper Belt's excited state (see Figure 
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4). However the viability of this scenario depends critically upon the timescales that govern 
the Belt's mass loss as well as the KBOs' velocity evolution. First note that accretion models 
show that as soon as the large R > 100 km KBOs form, further KBO growth is halted as they 
initiate an episode of more vigorous collisions that steadily grind the Belt's smaller members 
down to dust which radiation forces then transport away (Kenyon & Luu 1999; Kenyon & 
Bromley 2001). Accretion models show that the R ~ 100 km KBOs form at a ~ 35 AU 
over a r form ~ 1 x 10 7 year timescale, and that the Belt's subsequent collisional erosion 
occurs over a r crodc ~ 5 x 10 8 year timescale (Kenyon & Luu 1999; Kenyon & Bromley 
2001). And if Neptune's orbit had migrated substantially, the advancing 2:1 resonance 
would also have swept across the Main Belt, which further enhances the stirring as well the 
collisional/dynamical erosion. The resulting grinding and erosion thus makes it ever more 
difficult for the Belt to sustain waves as the initially massive Kuiper Belt mass is eroded by 
a factor of ~ 150, which also lowers the disk's \jq 10~ 3 . Note that the removal of the 
smaller KBOs also turns off the dynamical friction that was once keeping the particle disk 
quite thin. The gravitational stirring by the surviving larger KBOs is then free to raise their 
dispersion velocities c up to their surface escape velocity t> esc (Safronov 1972), which results 
in a thicker disk with \) ~ v csc /na ~ 0.02(i?/100 km) where R is the KBOs' characteristic 
size. Since this is substantially larger than the current Belt's P)q, it seems quite likely that 
this stirring shut off the waves while they were still of low amplitude (see Fig. 4) long before 
the Kuiper Belt was ground down to its present mass 5 . Consequently, apsidal and nodal 
waves were likely able to propagate throughout the Kuiper Belt during a timescale r wave that 
is bounded by the moment when the large KBOs first formed and when the Belt eroded 
away, i.e., Tf orm < r wavc < T erof j e - 



4.5. Comments on Studies of Spiral Waves in Galactic Disks 

It should also be noted that this implementation of the rings model does not account 
for any viscous damping of spiral waves which, as Hunter & Toomre (1969) point out, may 
be more important as bending waves approach a disk's outer edge. Unlike the sharp-edged 
disks employed here, a more realistic disk likely has a 'fuzzy' edge where the disk's surface 
density gently tapers to zero. Hunter & Toomre (1969) show that bending waves entering 



5 An exception to this assertion might occur if the reflection of nodal waves at an outer edge allowed the 
Belt to behave as a resonant cavity (Ward 2003). When the disk's mass and the wave's pattern speed are 
appropriately tuned, which can occur naturally as the Belt eroded and/or the as Neptune's precession rate 
varied due to its orbital migration, then a higher-amplitude standing wave pattern can result while the Belt 
persists in the tuned state. 
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such a lower-density zone tend to excite substantially larger inclinations there, and such sites 
will be considerably more susceptible to the viscous dissipation of these possibly nonlinear 
waves. 

Toomre (1983) suggests that disks having a tapered edge might also cause bending waves 
to stall there since the group velocity c g — > as the surface density a smoothly goes to zero. 
However this assertion was not confirmed by the rings model, which tapered a disk's outer 
surface density by multiplying it by the factor a/1 — {£ — Aa) 2 /£ 2 where Aa is the distance 
from the outer edge and £ is the tapering scale-length. These experiments adopt values of £ 
that are smaller, comparable, and larger than the bending wavelength Al, and in all cases 
the bending wave reflects at or near the disk edge, with considerably larger inclinations being 
excited in this tapered zone. 

5. Summary 

A model that rapidly computes the secular evolution of a gravitating disk-planet system 
is developed. The disk is treated as a nested set of rings, with the rings'/planets' time- 
evolution being governed by the Lagrange planetary equations. It is shown that the solution 
to the dynamical equations is a modified version of the classical Laplace-Lagrange solution 
for the secular evolution of the planets (Brouwer & Clemence 1961; Murray & Dermott 1999), 
with the modification being due to a ring's finite thickness h = c/n that is a consequence of 
the dispersion velocity c of that ring's constituent particles. Since the ring's finite thickness 
h softens its gravitational potential, this also softens the Laplace coefficients appearing in 
the Laplace-Lagrange solution over a scale h/a. 

It is shown that the Lagrange planetary equations admit spiral wave solutions when the 
tight-winding approximation is applied. There are two types of spiral density (or apsidal) 
waves, long waves of wavelength A^ = 27r//d|n/o;|a and short waves of wavelength As < 
2\/27r()a where f) = h/a is the disk's fractional thickness and uj is the angular rate at which 
the spiral pattern rotates. The simulations presented here show that the giant planets launch 
long waves at either a resonance in the disk or else at the disk's nearest edge, and that these 
waves propagate away until they reflect at the disk's far edge or else at a Q-barrier in 
the disk which resides where f) = f)g(a) where f)g(a) = 0.26fid\n/u\ is the maximum disk 
thickness that can sustain apsidal waves, with /id = naa 2 /M & being the normalized disk 
mass. Of course all of these findings may be derived from the stellar dispersion relation 
given in Toomre (1969) in the limit that the pattern speed u is much smaller than the 
disk's mean motion n. Nonetheless, it is satisfying to see that the theory of unforced apsidal 
waves is readily obtained from the Lagrange planetary equations; with a little more effort 
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the theory for forced apsidal waves [e.g., Ward & Halm (1998)] should also be recoverable. 

However new results are obtained for the nodal wave problem, which admits only a long- 
wavelength solution \ L to the planetary equations in the tight-winding limit. In particular, 
it is shown that these waves can stall, that is, the waves' group velocity plummets to zero as 
they approach a site in the disk where f) = 2.72[)q. If, however, these waves instead encounter 
a disk edge, they will reflect and return as long waves. In the limit that () — > 0, the results 
for nodal waves propagating in a infinitesimally thin disk is recovered (Ward & Hahn 2003), 
but note that the wave-stalling phenomenon does not appear in a f) = treatment of the 
disk. 

The rings model is also used to examine the propagation of apsidal and nodal waves that 
are launched by the giant planets into a variety of Kuiper Belts having a mass Mkb — 30 
M e (the estimated primordial mass) down to Mkb — 0.08 M e (which is ~ 40% of the 
Belts current mass estimate). In each simulation the giant planets deposit roughly the same 
fraction of their initial angular momentum deficits, ~ 0.5% and ~ 10% of the planets' L e 
and Li, into the disk in the form of spiral waves. And since the waves' angular momentum 
content is roughly the same in each simulation, the lower-mass Kuiper Belts thus experience 
higher-amplitude waves. Indeed, the waves seen in the Mkb < 0.2 M e simulations are of 
sufficient amplitude that they could in principle account for much of the dynamical excitation 
that is observed in the Kuiper Belt. However wave-action in a Mkb ~ 0.2 M e Belt also 
requires its fractional scale-height to quite thin, namely, f) < 10~ 3 . Most likely, apsidal and 
nodal waves were shut off, due to self-stirring by large KBOs as well as by other external 
perturbations, long before the Belt eroded down to its current mass, in which case the 
excitation by wave-action would have been quite modest. 

The rings model developed here has many other applications. One issue of great interest 
is to determine whether apsidal and nodal waves may be propagating in Saturn's rings. Of 
particular interest are the short apsidal waves since their detection could yield the ring's 
dispersion velocity c via a measurement of the short wavelength \s ~ 9c/ n. Although the 
ring particles' dispersion velocity is of fundamental importance to ring dynamics, it is less 
than well constrained at Saturn. Of course, the differential precession due to planetary 
oblateness also needs to be included in the model [c.f., Murray & Dermott (1999)] since 
this effect may actually defeat this form of wave-action. The rings model can also be used 
to examine the forced motions of a relatively massless but much thicker circumstellar dust 
disk like f3 Pictoris. The warps and brightness asymmetries seen in this system are usually 
attributed to secular perturbations exerted by an unseen planetary system, and the code 
developed here can be used to very rapidly explore the wide range of planetary parameters. 
This rings model will be used to study these and other problems in greater detail in the near 
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future. 
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A. Appendix A 

Differentiating the appearing in the / and g functions, equations (17), yields 

/ =2abf l2 + \a 2 (&$ - 6$) - 3^ 2 (2 + E 2 )b% (Ala) 

g =2(a 2 + !)(! + H 2 )b$ 2 - 3a (&$ + bf/ 2 ) + 3a: 2 tf 2 (2 + H 2 )~b^ 2 ~ \a 2 (&$ - 5$) 

(Alb) 

where if 2 = |(f) 2 + f)' 2 ). The following recursion relations, 

m&( m) = S « (Vfc* - VS+V) (A2a) 
(m + 1 - s)ati™ +1) =m(l + a 2 )(l + H 2 )b^ -(m + s- l)a&( m_1) , (A2b) 

can be used to simplify equations (Al) further. Equations (A2) are derived in Brouwer 
& Clemence (1961) for the case where H = 0. However the more general relations given 
above are readily obtained by replacing the combination 1 + a 2 appearing in the Brouwer 
& Clemence (1961) recursion relations with (1 + a 2 )(l + H 2 ). So for m — 1 and s = 3/2 
equation (A2a) is 

J& = §«(5&-5$) (A3) 

and 

5(2) _ l a (m _ m 

h/2 — 4" ^5/2 °5/2 / 

for m = 2 and s = 3/2 while equation (A2b) yields 



= -nib'M, - ;",) (A4) 



2(a 2 + 1)(1 + H 2 )~bW = ab^ + 3a6$ (A5) 

for m = 1 and s = 3/2. Inserting equations (A3-A5) into (Al) then yields 

f{a, f), fj') =0^/2 - 3a 2 # 2 (2 + tf 2 )^ (A6a) 

^(a, fj, fj') = - a6^ 2 + 3a 2 H 2 (2 + tf 2 )^. (A6b) 

B. Appendix B 

The symbolic mathematics software MAPLE has been used to write the needed softened 

Laplace coefficients b~i m \ equation (12), in term of complete elliptic integrals K and E. 
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Setting H 2 



±(h 2 + f)' 2 ), 7 = (1 + a 2 )(l + H 2 )/2a, and x = ^2/(7 + 1), then 



7TA/2a(7 + 1) 
_ 4[7K( X )-(7 + l)^(x)] 
71^2^(7 + 1) 

7T«( 7 - l) v /2a( 7 + 1) 

_ 2[-( 7 -l)K(x)+7^(x)] 

?ra(7 - l) v /2a(7 + 1) 

_ 2[-47(7-l)K( x ) + (47 2 - 

?ra(7 - l) v /2a(7 + 1) 

4[-(7-l)^(x)+4 7 ^(x)] 
"37r(2a) 5 /2( 7 + i)3/2( 7 _ 1)2 

4[- 7 ( 7 -i)ir( X ) + (7 2 + 3)E(x)] 
37T(2a) 5 / 2 ( 7 + l) 3 / 2 ( 7 - l) 2 



3)^(X)] 



where 



K{X) 



dt 



lo v/(l-t 2 )(l- X 2 t 2 
is the complete elliptic integral of the first kind and 



E(X) 



x 2 t 2 



1-t 2 



dt 



(Bla) 
(Bib) 
(Blc) 
(Bid) 
(Ble) 
(Blf) 
(Big) 

(B2) 



(B3) 



is the complete elliptic integral of the second kind. The series expansions for K(x) an d E(x) 
given in Abramowitz and Stegun (1970) then permit the rapid calculation of the softened 
Laplace coefficients without requiring a numerical integration of equation (12). However 
equations (Bl) can give unreliable results for extreme values of a due to numerical roundoff 
errors. In this case it is preferable to factor the (1 + a 2 )(l + H 2 ) = 2«7 term out of the 
integrand in (12) and expand the denominator for the case of large 7 ^> 1: 



~ h {rn) 



f m a n 
"(2«7) 



d(f> cos(m0) 



s s (g _i_ 1 ) 

1 + - COS0+ — — ^ COS 20+ ... 

7 27^ 



s+m 



(B4a) 
(B4b) 



where fo = 2, f\ = 2s, and / 2 = s(s + 1). Equation (B4b) usually gives the more reliable 
result for a <C 0.01 and for a 3> 100. 
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Another useful form for is obtained for regions where a = 1 + x where \x\ <C 1. In 
this 'local' approximation, the dominant contribution to the integral, equation (12), occurs 
where ^ < 1. Thus we can set cos0 ~ 1 — 2 /2 and a ~ 1 except where it appears 
as a — 1 = x, extend the upper integration limit to infinity, and set cos(m0) ~ 1 in the 
numerator [cf. Goldreich & Tremaine (1980)]: 



2/tt 



for s = 3/2 



:; J- 2// 2 

for s = 5/2 



(x 2 + 2# 2 )' 



C. Appendix C 

The time derivatives of L e , equation (26a), is 

^^mwXh^ + k^) (Cla) 

3 

= ^2^2 m j n2 j a ) A jk{hjk k - h k kj) (Clb) 

3 k^j 

= Yl Yl \ \ M 3 + k m ) ^K&^A ~ hkk^ (Clc) 

3 k^j & ' 3 / 

where g 3k = g(otjk, tyj, t)k)- Let dL e /dt — S± — S 2 where Si is the sum over the hjkk terms 
and S 2 the sum over the h k kj. Swap the j and k indices in S 2 so that it becomes 

k j¥^k 



Equation (18c) shows that g k j = OLj k g,j k , and with (rik/nj) 2 = a- k 3 {M Q +m fc )/(M +m,j), S 2 
becomes 

*-EEK£t^)"H*^ (C3) 

k jy^k x J 7 

which is Si since the sums obey V^ fc YE , fc = V^- X^j- Consequently dL e /dt = Si — S 2 = 0, 
and a similar analysis will also show that dLi/dt = 0. 
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wavenumber K 
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Fig. 1. — The upper figure is the dispersion relation oo e = Ke for apsidal density waves, 
and the lower figure in the dispersion relation Ui = e~ K — 1 for nodal bending waves. 
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Fig. 2. — A snapshot of the simulated Kuiper Belt's orbit elements (e, i, u>, f2) plotted versus 
semimajor axis a at time t = 2 x 10 6 years in a M KB = 10 M e Belt. The Belt's fractional 
half-thickness is f) = 0.0067. The dots along the left axes indicate Neptune's orbit elements. 
The Aa/a greyscale shows the disk's fractional surface density variations versus the polar 
coordinates (a, 0), estimated via equation (39). White/black zones indicate over/under dense 
regions where \Aa/a\ exceeds 0.63, while white/black zones in the other greyscale indicate 
the disk's latitude (5 — z/a above/below the invariable plane, with saturation occurring 
wherever \(3\ exceeds 0.86°. The online edition of this figure will be linked to an animated 
sequence of these figures that show the system's time-evolution. 



-38- 



precession period P, (years) 




eigenf requency g. (radians/year) 



10 7 10 6 10 



1 — i — i — i 1 1 1 1 1 1 — i — i — i 1 1 1 1 r 




eigenfrequency If, I (radians/year) 



-39- 



Fig. 3. — The numerous small dots show the individual elements in the disk-rings' eccen- 
tricity eigenvectors \Eji\, all plotted versus their eigenfrequencies <?j, as well as the disk-rings 
inclination eigenvector elements versus for the system shown in Fig. 2. The solid 
curves are \Eji\ and \Iji\ averaged over semimajor axes 45 < a,j < 55 AU. The upper horizon- 
tal axes are the corresponding precession periods 2n/gi and 27r/|/j|. The large dots indicate 
the eigenfrequency and eigenvector element that dominates the motion of each giant planet 
J=Jupiter, S=Saturn, U=Uranus, and N=Neptune. 
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Fig. 4. — The maximum eccentricities e max and inclinations i max versus semimajor axis a for 
simulations having a variety of Kuiper Belt masses M^b] see Section 4.2 for model details. 
Dots indicate the orbits of 340 KBOs observed over multiple oppositions, with orbits provided 
by the Minor Planet Center. The locations of the 3 : 2 and 2 : 1 resonances are indicated, 
and orbits above the grey curve are Neptune-crossing. 
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Fig. 5. — Maximum eccentricities and inclinations in a M^b = 10 M e disk having a fractional 
thickness fj = 0.0027,0.027,0.04,0.08, and 0.13. The black curves are the forced e's and i's 
that occur in a massless disk. The characteristic particles size associated corresponding 
to each disk thickness f) may be obtained by setting particle dispersion velocities equal 
to their surface escape velocity, which corresponds to particle radii of it! ~ 5000t) km ~ 
14, 140, 200, 400, and 650 km, respectively. 
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Fig. 6. — The upper plot gives the maximum eccentricities e max in a Mkb = 0.2 M e disk 
having thicknesses f} = 0.0015,0.0037,0.0074, and 0.011. The lower plot gives the average 
inclinations % that occur after the bending wave has stalled. The exception is the orange 
[) = 0.0015 curve; this wave does not stall, so this curve shows the maximum inclination 
i max - The black curves are the forced e's and i's that occur in a massless disk. 



